The jamming transition as probed by quasistatic shear flow 
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We study the rheology of amorphous packings of soft, frictionless particles close to jamming. 
Implementing a quasistatic simulation method we generate a well defined ensemble of states that 
directly samples the system at its yield-stress. A continuous jamming transition from a freely-flowing 
state to a yield stress situation takes place at a well defined packing fraction, where the scaling laws 
characteristic of isostatic solids are observed. We propose that long-range correlations observed 
below the transition are dominated by this isostatic point, while those that are observed above the 
transition are characteristic of dense, disordered elastic media. 
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A collection of spherical particles, interacting via a 
finite-range repulsive (contact) potential, unjams from 
a solid to a non-rigid state when being decompressed be- 
low a critical volume-fraction, (pc [H, This transition, 
which has been given the name "point J" , is accompanied 
by several interesting and nontrivial scaling relations in 
the solid phase [J, There, pressure and Hnear elastic 
shear modulus vanish as does the ratio of shear to bulk 
modulus. The average number of particle contacts jumps 
from a finite value zq at point J to zero just below the 
transition. The value of zq is compatible with Maxwell's 
estimate for the rigidity transition and signals the fact 
that at point J each particle has just enough contacts 
for a rigid /solid state to exist (which is called "isostatic" 
state) [J, |5|- Above the transition, additional contacts 
are generated according to the surprisingly universal law, 
Sz ^ Scj)^/'^. While the system moves away from its iso- 
static state the effective size of the remaining isostatic 
regions, I* ^ 5z~^ [1,01 has been argued to provide the 
diverging length-scale that endows point J with a certain 
"criticality" . Based on these findings, jamming has been 
regarded as a "mixed" transition that shares properties of 
both, discontinuous (jump in contact number) and con- 
tinuous phase- transitions (diverging length-scale) [1, Q . 

A different route to approach point J from the fluid 
phase has been used in the flow simulations of Refs. [H, 
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12j . Several scaling relations have been reported, 
some of which depend on model details. Some others 
seem to share the universality encountered in the solid 
phase. In contrast to the linear elastic properties in 
the solid phase, only little understanding about the flow 
properties and their relation to nearby point J has been 
achieved up to now. 

In this Letter, we employ a quasistatic simulation tech- 
nique that studies the borderline between fluid and solid 
state. As we will see, it combines both aspects, the elas- 
ticity of a solid and the flow of a fluid, in one simulation. 
This fact allows us to get insight into how point J, and 
its isostatic state, affects the flow properties close by. In 
particular we show that jamming (i.e. the development 



of a finite shear stress) as probed by shear flow should 
be viewed as a continuous transition, with isostatic ef- 
fects showing up in the flow properties primarily below 
the transition (0 < 0c). In contrast, at volume- fractions 
above 0c the flow is characteristic for amorphous, but 
well connected, materials. 

The quasistatic simulation probes the flow of the sys- 
tem in the limit of vanishing shear rate, 7 — > (see 
Fig. [1]). On lowering the volume-fractions from values 
above (pc point J is therefore approached along the yield- 
stress line. The resulting "yield-stress-flow" is generated 
by a succession of equilibrated solid states that, as we 
will see, carry the signature of the nearby isostatic state. 
At point J the yield stress vanishes such that at lower 
volume- fractions, < 0c, the simulation follows the 0- 
axis and the system flows at zero stress. This corresponds 
to the limiting case, 7 0, of normal fluid flow in the 
Newtonian regime [10[. 



O 



FLUID 



yield-stress 

SOLID 




FIG. 1: Trajectory taken by quasistatic simulations in state- 
space of shear stress cr and volume-fraction 0. (Dotted lines) 
contours of constant strain rate 7. The simulation corre- 
sponds to 7 ^ and thus follows the yield-stress line, ay{<j)). 
Previous approaches either probe linear elasticity of the solid 
(decompression at cr = 0) or the steady-state flow of the fluid. 

Our system consists of soft spherical particles that 
interact with an harmonic contact interaction with spring 
constant k. As in Refs. [ll,[ll,[l3| we use a 50 : 50 mixture 
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of particles with radius a and 1.4a in two dimensions. 
In the simulations the volume fraction cf) is controlled, 
while stresses are allowed to fluctuate. Primary output 
are the pressure, p, and (shear-)stress ct as a function of 
imposed strain, 7 (see Fig. [5]). To implement the shear, 
variable Lee-Edwards boundary conditions are used. The 
simulation proceeds by minimizing the total potential en- 
ergy (using conjugate gradient techniques [ij]) after each 
affine change in the boundary conditions and particle co- 
ordinates, (A7 = 5-10~^). Thus, as the energy landscape 
evolves under shear the system always remains at a local 
energy minimum, with all forces fully equilibrated. 

As can be seen from Fig. [2] a typical feature of qua- 
sistatic stress-strain relations is the presence of clastic 
branches, where the stress grows linearly with strain. 
This reversible clastic loading is terminated by irre- 
versible plastic events during which the stress drops 
rapidly and energy is dissipated. The succession of elas- 
tic and plastic events defines the flow of the material just 
above its yield-stress cr-y ((/>). On lowering the volume- 
fraction the average stress- level decreases until, below (j)c, 
it vanishes and the system flows at zero stress. For in- 
termediate volume-fractions, close to (j}c, one infers from 
Fig. [21 that the stress-signal is highly intermittent show- 
ing a coexistence between the two states of yield-stress 
and zero-stress Newtonian flow. Note, the striking simi- 
larity with the stress-strain curves of the experiments of 
Behringer et al. (see e.g. fisf). 



"liquid-like" to "solid-like". 




FIG. 2: Stress-strain relation for a sample of 2500 particles 
at two different volume fractions, (j) — 0.8470, 0.8433. 

We have counted the number of jamming-events that 
lead from the zero-stress flow to the yield-stress flow 
at finite stress, and back (see Fig. [3]). At high volume 
fractions the system always flows at finite stress so no 
event occurs. The same is true at low volume fractions 
where the system always flows at zero stress. Thus, we 
expect a maximum at an intermediate volume fraction 
0c where the system is highly unstable and makes fre- 
quent transitions between the two types of flows. For 
increasing system size wc find the width of the coexis- 
tence region to decrease, presumably vanishing in the 
thermodynamic limit. The extrapolated peak position, 
0c(oo) = 0.8433 ± 0.0003, can thus be used to define 
a critical density. At this volume-fraction the system 
develops a yield-stress and the flow changes from being 




FIG. 3: Number of jamming events (defined as states where 
stress first exceeds a threshold value, athrosh = 2e — 07) as 
function of volume fraction (j> ^nd system size A^. Inset; 
Finite-size scaling indicates a vanishing width and a diverging 
height as A'' 00. 

The results obtained here are similar to the finite-size 
scaling in Ref. [l[. There is an important difference, how- 
ever. In that study statistical information is gathered by 
generating initial states that are as random as possible. 
The observable is the probability that such a randomly 
chosen state is jammed. In contrast, here we are con- 
cerned with states that are connected by the trajectory 
of the system itself. We are quantifying a dynamical tran- 
sition rate at which the system, in the course of shear, 
undergoes a jamming transition. 

One may view the equilibrated finite-stress states vis- 
ited during the flow (0 > 0c) as an ensemble of jammed 
states that is generated, not by an external protocol, but 
by the energy landscape itself. In this spirit we have an- 
alyzed the elastic branches and found good agreement 
with previously identified anomalous scaling properties 
characteristic of marginally rigid solids. In particular, 
we found the shear modulus g to scale with the square 
root of pressure, g ~ p^^^ (not shown). Note, however, 
that this shear modulus is not the linear elastic modulus 
of the system g{a = 0) = gun, as discussed in Ref. [l[. 
Rather, it is defined at (or just below) the yield-stress, 
g = g{(7 = (7y). As the yield-stress itself is pressure de- 
pendent (wc find Gy ~ p) the modulus may in fact show 
a complicated dependence on pressure, p. This is not the 
case, however, if we assume g to take the scaling form, 
g[(J,p) = p^/'^F[a/ay{p)). At the yield stress, the argu- 
ment of the scaling function is a constant and we recover 
g ~ p^/'^ just as at zero stress. Note, that the presence of 
a scaling function F is to be expected close to point J. It 
has furthermore explicitly been evidenced in the elastic 
network model of Ref. [16| . 

Another hallmark of marginally rigid solids is the scal- 
ing of the contact number, z. In agreement with the sec- 
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nario under decompression we find, = z — zq ~ y^^"^ 
(see Fig. |4]). The number of contacts at zero pressure, 
zq w 3.8, is slightly smaller than the isostatic value, 
as rattlers have not been accounted for [2^. The scal- 
ing has been shown to result from a competition be- 
tween two terms in an expansion of the elastic energy. 



W ~ — pu^ [l7|. In this expression, U|| denotes the 



change in compression of a particle contact, while u± 
relates to the rotation of the two particles around each 
other. It has been shown (@, 0]) that both displacement 
components are related by ity ~ u±Sz, such that the en- 
ergy can be written as, W/u'^ kSz"^ — p. Thus, both 
terms are of the same size, when Sz'^ ~ p, which is the 
scaling relation observed numerically. 
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FIG. 4: Contact number z as function of pressure (left inset) 
and volume- fraction (j> (main panel). Symbols arc as in Fig. [3] 



Right inset: scaling of gap width, Sz 



; 3.8— z, with system- 



above and below suggest that different mechanisms are 
responsible for the formation of contacts. Indeed, the en- 
ergetic competition that controls the number of contacts 
in the jammed system is absent below (/)c, where con- 
tact formation is of purely geometric origin (the excluded 
volume of the particles). Interestingly, this also implies 
that the isostatic length-scale I* ~ scales differ- 

ently with volume- fraction above and below jamming. 
To define I* below (pc, note, that an isostatic cluster can 
be equilibrated by blocking its " surface" degrees of free- 
dom. The scale I* is then defined (similarly as above 0c 
SH) size of clusters, in which the number of extra 

degrees of freedom in the bulk (due to missing contacts, 
5zl*'^) is comparable to the number of surface degrees of 
freedom, 

To summarize this part, we find the ensemble of solid 
states at the yield-stress, cry(0), to display the scahng 
laws characteristic of marginally rigid solids. As the 
contact number does not show a discontinuity at (j)c we 
should view jamming under shear as a conventional con- 
tinuous transition, in contrast with the mixed character 
ascribed to the transition seen under decompression. 
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: 0.844, 0.843. 



In decompression simulations, the number of contacts 
displays a discontinuity at c/)c and jumps from the iso- 
static value, zq, to zero. This reflects the absence of 
any dynamics that would rearrange the particles once 
they have lost contact. In contrast, steady shear as 
studied with the quasistatic simulation leads to struc- 
tural rearrangements and thus to particle contacts even 
at low volume fractions below (pc- As can be seen in 
Fig. m there is still a gap in the number of contacts 
between jammed finite-stress (upper branch) and zero- 
stress (lower branch) configurations. However, this is just 
a sampling artifact and reflects the tendency for jammed 
configurations to have z > zq, while zero-stress states 
generally have z < zq. As can be verified in the inset, 
the width of the gap is reduced with increasing the sys- 
tem size and seems to vanish in the thermodynamic limit, 
where only states at zq are sampled. 

Furthermore, we find (see Fig. [4]) that below (pc 
the number of contacts simply increases linearly with 
volume-fraction, 5z ~ 5(p, in striking contrast with the 
square-root behavior found above (f>c- The resulting cusp 
at 0c is illustrated in the figure. The different scaling 



FIG. 5: Non-afRne displacements A (A'' = 900) taken from: 
average over all configurations (red plus), restricted aver- 
age over fiowing states (circles), jammed states (square) and 
jammed states without plastic events (triangles). (Small 
closed symbols) flowing/ jammed states at A'' = 1600. 

This continuity is also present on the single-particle 
level as we will discuss now. For non-interacting particles, 
the particle displacements are affine, u = 'yxey (for shear 
in y-direction) . Interactions lead to additional non-affine 
motion, in particular to a non-zero x-component, Wna = 
u ■ ex. In Fig. [5] we display A = {u^g}^^"^ as a function of 
volume fraction, 0. 

We find perfect continuity of the particle displacements 
across (pc only if we exclude the plastic events from the 
average and consider the displacements either in the re- 
versible elastic states (triangles) or the states at zero- 
stress (circles). Plastic events are, in general, violent re- 
arrangements that span the whole system. Consequently, 
these events dominate the amplitude of the mean-square 
displacement in the finite-stress states (squares). 

From the continuity of the displacements we conclude 
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that the mechanism that leads to nonafhne motion in the 
ensemble of solid states above 4)c is the same as that in the 
fluid flow just below (f)c- The difference being that above 
4>c particle motion leads to the build-up of stress and 
its subsequent release in the plastic events. Apparently, 
this does not affect the magnitude of the (elastic) non- 
afhne displacements. However, it does show up in their 
correlations. 
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at small strain-rates: the flow properties at small (and 
zero) stress and volume-fractions below <f>c are governed 
by the vicinity to the isostatic state. This is evidenced 
by a correlation length that diverges on approaching the 
isostatic point J from below. As a consequence, one may 
view the flow as due to the rearrangements of a liquid 
of marginally rigid, isostatic clusters. Close to point 
J a cross-over takes the system to a different regime 
{4> > 4>c)i where the flow properties are dominated by 
long-range elastic interactions. This regime, thus, reflects 
the physics of amorphous, but well connected, materials, 
where flow is due to the irreversible rearrangements of 
liquid-like defects in a solid matrix. 

The authors acknowledge fruitful discussions with 
Pinaki Chauduri and Ludovic Bcrthier, as well as thank 
the von-Humboldt Feodor-Lynen, the Marie-Curie Eu- 
rosim and the ANR Syscom program for financial sup- 
port. 



FIG. 6: Correlation function G{y) {N = 2500) for various 
volume-fractions above and below ■ For ^ correlations 
do not decay to zero, probably due to the peridiodicity of 
the simulation box. Inset: Length-scale ^ as defined by the 
position of the minimum. Dashed lines are ~ Scp^^'^ and 
~ 50-^ ((^c = 0.8433). 

Following Ref. [l^ we study the correlation function 
G(r) = (Mna(r)una(O)) takcu at r = e^y thus point- 
ing along the shear direction. As in Ref. [10| we find a 
minimum of the correlation function in the freely flowing 
phase {(j) < 4>c)- The corresponding length-scale ^{Scf)) 
grows with approaching the critical volume fraction from 
below. Unfortunately, extracting a reliable value for the 
correlation length exponent v is hampered by the re- 
stricted range of only one order of magnitude. Values 
may range from i^ = 0.8...I.0as indicated in the inset of 
Fig. [6] (no fit). Note, that a data fit would be highly sen- 
sitive to the value of (j>c used. This may also be the rea- 
son why previous studies report smaller values, 0.6 ... 0.7 
[lol . [lit ■ With an exponent of i/ = 1 this may indeed be 
the isostatic length-scale, ^ = P ^ Sz^^ ^ S(j)~^, thus, 
indicating the presence of isostatic clusters that grow on 
approaching (j)c from below. 

A quite different behavior is found in the regime above 
0c- Here G{y) does not have a minimum and always 
decays monotonously. The relevant length-scale is the 
system-size, ^ ^ L. This is consistent with findings 
in Ref. [l3| where the lack of observable length-scale is 
explained with the dominance of long-range elastic cou- 
plings that place yield-stress flows automatically into a 
system-size dominated regime. 

These results suggest the following picture for the flow 



[1 
[2: 

[6: 

[7 

[8 

[9: 
[10: 

[11 

[12 

[13 

[14 
[15 

[16 

[17 

[18 

[19 
[20 



C. S. O'Hern, L. E. Silbert, A. J. Liu, and S. R. Nagel, 
Phys. Rev. E 68, 11306 (2003). 

H. A. Makse, N. Gland, D. L. Johnson, and L. M. 

Schwartz, Phys. Rev. Lett. 83, 5070 (1999). 

T. S. Majmudar, M. Sperl, S. Ending, and R. P. 

Behringer, Phys. Rev. Lett. 98, 058001 (2007). 

C. F. Moukarzel, Phys. Rev. Lett. 81, 1634 (1998). 

A. V. Tkachenko and T. A. Witten, Phys. Rev. E 60, 

687 (1999). 

M. Wyart, S. R. Nagel, and T. A. Witton, Europhys. 
Lett. 72, 486 (2005). 

M. Wyart, L. E. Silbert, S. R. Nagel, and T. A. Witten, 
Phys. Rev. E 72, 51306 (2005). 

S. Henkes and B. Chakraborty, Phys. Rev. Lett. 95, 
198002 (2005). 

J. M. Schwarz, A. J. Liu, and L. Chayes, Europhys. Lett. 
73, 560 (2006). 

p. Olsson and S. Teitel, Phys. Rev. Lett. 99, 178001 

(2007) . 

T. Hatano, J. Phys. Soc. Jpn. 77, 123002 (2008). 

P.-E. Peyneau and J.-N. Roux, Phys. Rev. E 78, 011307 

(2008) . 

T. K. Haxton and A. J. Liu, Phys. Rev. Lett. 99, 195701 

(2007). 

|http: / /lammps.sandia.gov/index.html[ 
R. P. Behringer, D. Bi, B. Chakraborty, S. Henkes, and 
R. R. Hartley, Phys. Rev. Lett. 101, 268301 (2008). 
M. Wyart, H. Liang, A. Kabla, and L. Mahadevan, Phys. 
Rev. Lett. 101, 215501 (2008). 

W. G. EUenbroek, E. Somfai, M. van Hecke, and W. van 
Saarloos, Phys. Rev. Lett. 97, 258001 (2006). 
J. A. Drocco, M. B. Hastings, C. J. Reichhardt, and 
C. Reichhardt, Phys. Rev. Lett. 95, 088001 (2005). 
C. E. Maloney, Phys. Rev. Lett. 97, 035503 (2006). 
We have also performed simulations with rattlers defined 
as particles with zero contacts. The resulting estimate for 
zo is much closer to the expected value of 4. 



